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Abstract 

The bubbles involved in sonochemistry and other applications of cavitation oscillate inertially. A correct estimation of 
the wave attenuation in such bubbly media requires a realistic estimation of the power dissipated by the oscillation of 
each bubble, by thermal diffusion in the gas and viscous friction in the liquid. Both quantities and calculated numerically 
for a single inertial bubble driven at 20 kHz, and are found to be several orders of magnitude larger than the linear 
prediction. Viscous dissipation is found to be the predominant cause of energy loss for bubbles small enough. Then, 
the classical nonlinear Caflish equations describing the propagation of acoustic waves in a bubbly liquid are recast and 
simplified conveniently. The main harmonic part of the sound field is found to fulfill a nonlinear Helmholtz equation, 
■where the imaginary part of the squared wave number is directly correlated with the energy lost by a single bubble. For 
low acoustic driving, linear theory is recovered, but for larger drivings, namely above the Blake threshold, the attenuation 
coefficient is found to be more than 3 orders of magnitude larger then the linear prediction. A huge attenuation of the 
wave is thus expected in regions where inertial bubbles are present, which is confirmed by numerical simulations of the 
'nonlinear Helmholtz equation in a f D standing wave configuration. The expected strong attenuation is not only observed 
but furthermore, the examination of the phase between the pressure field and its gradient clearly demonstrates that a 
traveling wave appears in the medium. 

fCeywords: Acoustic cavitation. Bubble dynamics. Propagation in bubbly liquids. Wave attenuation 
PACS: 43.25.Yw, 43.35.Ei, 43.25. Gf 
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1. Introduction 

The complexity and large variety of spatial and tem- 
]poral scales involved in acoustic cavitation make difficult 
■the derivation of a full theoretical model, accounting for 
the coupled effects between the bubble field and the sound 
•field. Nevertheless, considerable progress has been made 
in the last decade. Theoretical studies in the context of 
single bubble sonoluminescence have allowed to restrict 
the ambient size of the bubbles in the micron range, owing 
ito surface instabilities [H-Q- This has been confirmed by 
20 kHz experiments, both in single bubble fs'] and multi- 
bubble configurations (see Ref. [9j] and references therein). 

Besides, a large collection of experimental observations 
■have revealed that radially oscillating bubbles in high- 
intensity acoustic fields tend to self-organize into bubble 
structures, which shapes depend on the experimental con- 
figuration, with possibly two structures or more appearing 
simultaneously in different zones of the liquid [ill|. The 
shape of such structures is strongly correlated with the 
fundamental issue of the translational motion of the bub- 
bles under Bjerknes forces, which have been reconsidered 
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in the context of strongly nonlinear inertial radial oscilla- 
tions and traveling waves 

On this basis, the action of the acoustic field on the 
organization of inertial bubbles has been satisfactorily de- 
scribed in various configurations by particle models [7, f ^ 
by calculating the forces exerted on the bubbles di- 
rectly from their nonlinear dynamics. Assuming a simple 
shape of the sound field, some bubble structures have been 
remarkably caught by this method. However, the correct 
prediction of other structures was found to be more diffi- 
cult, mainly because, as suggested by Mettin [9], the local 
sound field might have a complicated shape, which cannot 
be inferred without describing correctly the acoustic field 
in the medium. 

The backward effect of inertial bubbles on the propa- 
gation of acoustic waves remains mainly unexplored. The 
main physical effects of the bubble radial oscillations on 
sound waves can be easily understood qualitatively. Bub- 
bles are mechanical oscillators so that wave dispersion is 
expected. They oscillate non linearly for large amplitude 
drivings, so that waves should be nonlinear. Finally, they 
dissipate mechanical energy by various processes, which 
should produce wave attenuation. The problem has been 
attacked in the early work of Foldy [l7| who considered 
linear scattering of waves by an arbitrary statistical distri- 
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bution of scatterers, and obtained a linear dispersion re- 
lation. The application of this theory to the specific case 
of linear sound waves in bubbly liquids has been consid- 
ered in Refs. [H, [l^. A key feature in Foldy's approach 
is that for a sufficiently dilute bubbly mixture, each bub- 
ble behaves as if it were excited by the statistical average 
pressure field, which allows to cast aside the difhcult issue 
of bubbles pairwise interaction. An intuitive justification 
of this approach can be found in Refs. 2^, 21 1. 

The assumption of small amplitude waves has been re- 
laxed by lordansky 22| and simultaneously by van Wijn- 
gaarden psi . [23 ] by a semi-empirical volume-averaging of 
the bubbly liquid equations, which is closed by a Rayleigh 
equation, where, as suggested by Foldy's work, the driv- 
ing pressure term is the local average pressure field. The 
model obtained has allowed the study of nonlinear dis- 
persive waves. The latter are classically described by the 
Korteweg-de Vries equation (isj , and the reduction of van Wi- 
jngaarden model to the latter for moderate amplitudes has 
been studied by various authors both theoretically [24 . 
[2^ and experimentally 29l-3lj|. 

The popular Caflish model (32| is a rigorous general- 
ization of Foldy's theory to the nonlinear case and yields 
a simplified version of van Wijngaarden model, as far as 
the bubbly liquid is dilute enough. The latter hypothesis 
has the important corollary that the mean velocity of the 
mixture is infinitely small, so that the momentum conser- 
vation equation coincide with the one of linear acoustics 
[see Eq. ^ below]. A physical discussion of the latter 
feature can be found in Refs. [20l. [2l|. 

Under the linear approximation, the Cafiish model re- 
duces to the famous dispersion relation of Foldy, which 
can be extended to calculate a wave attenuation coeffi- 
cient, accounting for dissipation by a linearly oscillating 
bubble [2^ and to polydisperse bubbles size distributions. 
Linearization allows a simple description of the sound field 
by an Helmholtz equation, and has been used in studies 
of the coupling between wave propagation and the bubble 
field. The gain obtained by simplifying the wave equation 
allows a complex description of its coupling with the bub- 
ble population evolution, spatially and along the size axis. 
Following such an approach, Kobelev & Ostrovski 33] have 
proposed an elegant model of self-action of low amplitude 
sound waves in bubbly liquids, accounting for the bub- 
ble drift under the action of primary Bjerknes forces and 
bubble coalescence favored by secondary Bjerknes forces. 
Although the wave equation in this study was linear, the 
global model was nonlinear, owing to the dependence of 
the wave number on the varying bubble density, which 
conversely evolves non linearly with the sound field. Spe- 
cific solutions under different hypothesis could catch the 
experimentally observed self-transparency, self-focusing of 
sound waves in bubbly liquids, and destabilization of ho- 
mogeneous bubble distributions. The latter instability has 
also been demonstrated in Ref. 34] by a similar approach, 
but involving a slightly different physics. 

The attenuation of sound waves by oscillating bubbles 



remains normally weak for linear waves, except when the 
bubbles are close to the resonant size [il0,|3i|, which IS 
the main cause of sound extinction considered in Ref. (ssj . 
Since low-frequency inertial cavitation involves bubbles 
much smaller than the resonant size the use of the 
linear theory of Ref. predicts an abnormally low at- 
tenuation, compared to experimental data (36| . This is 
not astonishing since inertial bubbles typically suffer a ten- 
fold expansion of their radius and are expected to dissipate 
more energy than predicted by linear theory. Despite the 
latter restriction, the linear dispersion relation has often 
been used to predict attenuation of strong cavitation fields, 
because it allows the description of the problem by a linear 
Helmholtz equation, which is easy to solve, and allows har- 
monic response simulations 37h40|. Moreover, the use of 
a complex wave number provided by the linear dispersion 
relation in an Helmholtz equation somewhat masks the 
fact that one physical origin of wave attenuation by the 
bubbles is the energy dissipated by the latter. The latter 
point has been nicely addressed by Rozenberg [4l|, who 
restated the problem of attenuation of a traveling wave by 
a cavitation zone in terms of energy conservation, without 
resorting to the linear hypothesis. The latter study made 
use of an empirical expression, fitted on experimental re- 
sults, between the power dissipated by cavitation bubbles 
and the wave intensity. Doing so, realistic attenuated in- 
tensity profiles near the emitter could be calculated simply, 
and experimentally observed self-attenuation of the wave 
could be accounted for. 

The last remarks suggests that the relaxation of the 
linear hypothesis is necessary to correctly predict attenu- 
ation by inertial bubbles, so that one should revert to the 
original fully nonlinear form of the Cafiish model. How- 
ever, although valid for any wave amplitude, the latter re- 
mains intractable for large multi-dimensional geometries, 
since it requires time-dependent simulations, and presents 
convergence problems in the range of inertial cavitation. 



even in ID [42|, |43|. Thus, an intermediate model, simple 
enough to be numerically tractable, but properly account- 
ing for the true energy dissipation by inertial bubbles, is 
necessary. 

The motivation of this work is the derivation of such a 
reduced model, and can be viewed as a systematic formal- 
ization of Rozenberg's approach 4l|, based on the non- 



linear Caflish model. The present paper extends the ideas 
formerly presented in Ref. (44] and is organized as follows. 
In section 2, we recast the fully nonlinear Caflish equations 
into a mechanical energy balance equation, where we ex- 
press explicitly the energy lost by the bubbly liquid on 
average over an oscillation period, as functions of period- 
averaged quantities of a single bubble dynamics. This en- 
ergy loss is then computed numerically, by simulating a 
bubble radial dynamics equation over a typical parameter 
range, including the range of inertial cavitation involved 
in cavitation and sonochemistry experiments. In section 3, 
we then seek a reduction of the Caflish equations for the 
main harmonic component of the acoustic fleld, involving 
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the energy dissipation calculated in section 2. Finally, in 
section 21 the resulting nonlinear Helmholtz equation is 
solved numerically in a ID configuration, and a detailed 
analysis of the obtained wave profiles is performed. The 
implications of the present results on the primary Bjerknes 
forces and 2D simulations of classical experimental config- 
urations are deferred in a companion paper. 

2. Theory 

2.1. Caflish equations 

The Caflish model [32| describes the propagation of an 
acoustic wave of arbitrary amplitude in a bubbly liquid 
described as a continuum, which means that the radial 
oscillations of all the bubbles pertaining to an elementary 
small volume of mixture located at a spatial point r can be 
described by a continuous spatio-temporal radius function 
R{r,t). The first two equations of the model correspond 
to mass and momentum conservation in the mixture: 
I dp dp 
+ ^-" = W 



PlCi 



Pi 



(1) 



9v 

dt 



Vp = 0. 



(2) 



In the above equations, p{r, t) is the acoustic pressure field, 
v(r, t) the velocity field, pi the liquid density, q the sound 
speed in the liquid, and /3{r,t) is the instantaneous void- 
fraction, which, assuming a mono-disperse distribution of 
the bubbles, can be defined by: 



/3(r,i) =7V(r)-^i?(r,i)', 



(3) 



where A^(r) is the local bubble density. The latter is as- 
sumed time-independent, or at least almost constant on 
the time scale of the oscillations. Despite the set of equa- 
tions (dJ and ^ is very similar to the equations of lin- 
ear acoustics, the presence of the right-hand-side term of 
Eq. ([Ij renders the whole model nonlinear. Following the 
procedure classically used for linear acoustics, these two 
equations can be easily recast into an equation of energy 
conservation, by multiplying by p and (U) by v: 



d_ 

dt 



1 

2picf 



-.Piv 



V. (pv) = Np 



dV_ 

dt ' 



(4) 



where V{y, t) denotes the instantaneous volume of the bub- 
bles located at r. The time derivative in the left-hand-side 
(LHS) of this equation represent the time- variations of the 
acoustic energy density, which is the sum of kinetic energy 
and potential compressional energy of the pure liquid. The 
second LHS-term is the divergence of the acoustic intensity 
pv. The right-hand-side (RHS), which would be zero for a 
linear wave propagating in the pure liquid, represents the 
mechanical power exchanged between the acoustic wave 
and the bubbles. As will be seen below, part of this en- 
ergy is irreversibly dissipated along the radial oscillations 
of the bubbles, which is the physical origin of the acoustic 
wave attenuation. 



2.2. Bubble dynamics 

The bubble radial motion equation can be described by 
a radial dynamics equation. The Cafiish model in its orig- 
inal form uses a inviscid Rayleigh-Plesset equation with 
isothermal behavior of the bubble, in which the infinite 
driving pressure field is the mean local acoustic pressure 
field p(r, t). In the present study, we want to examine the 
energy dissipation by heat transfer between the bubble 
interior and the liquid, and by viscous friction in the ra- 
dial motion of the liquid around the bubble. We therefore 
leave the bubble pressure Pg unspecified for now, and add 
the classical viscous term in the Rayleigh-Plesset equation. 
Besides, since surface tension plays a preponderant role in 
inertial cavitation (isl - lisj , we also added the correction ac- 
counting for the latter effect, so that the bubble dynamics 
is given by: 



pi[RR+ -W = 



where fii is the liquid dynamic viscosity, and a the surface 
tension. All the quantities i?, Pg and p in this equation 
are spatio-temporal fields, depending on both r and so 
that the time derivatives represented by over-dots in this 
equation must be understood as partial derivatives d/dt 
at r constant. We did not add any corrections accounting 
for liquid compressibility, in order to keep a reasonably 
simple model. We defer the discussion of this choice to 
the conclusion section. 

For further use in the paper, we recall that when a bub- 
ble is driven by a sinusoidal pressure field p = po[\ — P* sin(27r/i)] 
around the ambient pressure po, its oscillations become 
inertia-controlled and involve a strong collapse when the 
dri ving pressure amplitude is above the Blake threshold 
ii-IiTi: 



1 



271 + 5 



1/2 



(6) 



where S = 2(t/(po^o) is the dimensionless Laplace ten- 
sion and Rq the bubble ambient radius. Such an oscilla- 
tion regime, historically termed as "transient cavitation" , 
is now classically named as "inertial cavitation" . 

2.3. Energy dissipation per bubble 

In order to get an energetic interpretation of the bubble 
radial motion, equation ([5]) can be multiplied by the time 
derivative of the bubble volume dV/dt, and noting that 
pi{RR + \R^) X dVjdt is the time-derivative of the radial 
kinetic energy of the liquid Ki = 2ttpiR'^R^, we obtain: 



d_ 

dt 



[Ki + A-KR^a) 



-WTTpiRlf 



dV 



dV 



(7) 



This equation is strictly equivalent to the Rayleigh equa- 
tion, and is the expression of the theorem of kinetic energy 
applied to the liquid surrounding the bubble. The paren- 
thesis in the LHS of ([7]) represents the sum of the kinetic 
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energy of the radially moving liquid and the interfacial 
potential energy. 

The first term in the RHS of Eq. ([7]) is the power irre- 
versibly lost by internal viscous friction within the liquid 
as it moves radially. 

The second term in the RHS is the power transferred 
from the acoustic field to the liquid surrounding the bub- 
ble, and can be viewed as the energy source available to 
drive the bubble oscillations and the radial motion of the 
liquid around. When multiplied by the number of bubbles 
per unit-volume, this term is similar to the right- hand- side 
of Eq. (|4]) with the opposite sign, which clearly indicates 
how energy is transferred between the driving acoustic field 
and the radially oscillating bubble. 

Finally, the last term in the RHS of (0 is the mechan- 
ical power done by the gas on the liquid, and could be 
expressed as the time-derivative of a compressional energy 
—dEp/dt in the case of a barotropic relation between the 
bubble pressure Pg and volume V (for example assuming 
an isothermal [32| or adiabatic evolution of the gas) . How- 
ever, in the general case where heat flows irreversibly be- 
tween the bubble interior and the liquid, this term cannot 
be expressed as the time-derivative of a potential function, 
and we now detail how this term is linked to dissipation of 
energy over a whole oscillation cycle of the bubble. 

In what follows, we will assume periodic oscillations of 
all the fields. Averaging Eq. (O over one cycle, the time- 
derivative in the left side cancels and we get: 



dV_ 
dt 



= n 



th 



Hv 



(8) 



where the two bubble dynamics-dependent average quan- 
tities Hth and Hv read: 



H 



th 



1 dV , 



H„ 



1 

r7o 



lenfiiRR^ dt. 



(9) 
(10) 



The quantity Hv defined by pil| is clearly positive, and is 
the period-averaged power loss by viscous friction in the 
liquid. 

A clear interpretation of Hth can be obtained by ap- 
plying the first principle of thermodynamics to the whole 
bubble content, which yields integral ^ as: 



H 



th 



TJo dt 



dt 



1 



Q dt, 



(11) 



where Ug and Kg depict the internal energy and kinetic 
energy, respectively, of the whole gas in the bubble, and Q 
is the heat gained by the bubble over one cycle. The first 
integral in the right-hand-side of pT|) is zero for a periodic 

motion, so that Hth = — (q^ is just the net heat lost by 
the bubble over one oscillation cycle. 

Equation ([5]) has therefore the following physical mean- 
ing: the energy transferred by the acoustic field to the 



bubble over one acoustic period is dissipated by two pro- 
cesses: the heat flow from the bubble toward the liquid 
and the viscous friction in the liquid radial motion. 

The integrals Hv and Hth can be evaluated numeri- 
cally by solving the bubble dynamics equation ([S]) for an 
arbitrary single bubble of ambient radius Rq excited by a 
sinusoidal forcing p — pq[1 — P* sin(27r/t)], possibly vary- 
ing the acoustic parameters P* and /, the bubble ambient 
radius Rq, and the properties of the liquid and the gas. In 
this paper we will restrict to air bubbles in water at ambi- 
ent pressure excited at 20 kHz and take: po = 101300 Pa, 
pi = 1000 kg/m3, pi = 10-3 Pa.s, a = 0.0725 N.m"!. The 
bubble ambient radius Rq and driving pressure amplitude 
P will be varied within a range of interest. More results 
involving, among others, the effect of the frequency and 
the type of gas will be given elsewhere [HO] • 

Since Hth represents the net heat flow leaving the bub- 
ble, thermal diffusion in the bubble interior must be prop- 
erly accounted for in our simulations, at least in an ap- 
proximate manner. To that aim, the bubble interior is 
modeled by an approximate energy conservation equation 
based on a thermal diffusion layer 51 -H^]. Water evap- 
oration and condensation at the bubble interface is also 
taken into account by a similar method, as described in 
the same references. The latter refinement may be im- 
portant since the presence of water vapor in the bubble is 
known to decrease the temperature collapse [F?! and there- 
fore influences the estimation of Q. To solve the bubble 
radial dynamics, the variables are non-dimensionalized by: 



t* 



Lot, 



R* 



R 



Pa 



Pg_ 

Pa' 



H* 



H. 



,th 



,th 



PqVouj' 



and the dimensionless dissipation functions H^j^ and H* 
are calculated numerically by: 



H 



th 



n* 



6 LOpi 

Ti" Pa Ja 



2a 

PaRa 



R* 



dR* 



' dV* 

Pr,—, dt 

dt* 



dt* 



(12) 



(13) 



For comparison purposes, we recall that assuming linear 
oscillations of the bubble, the equation of radial dynamics 
can be linearized by setting 



R*{t) = l + ^{Xe 



it' 



c.c. 



(14) 



where the complex amplitude can be obtained analytically, 
accounting rigorously for thermal effects 28, 5^, E^l- In- 
troducing (HH) in ([15]) -(dH), H^i^ and H* can be obtained 
analytically, and we obtain: 



h: 



th.hn 

TT* 



1 



Po 



2a 

PaRa 



S($,)|^| 



X\ 



(15) 
(16) 



where <l>g is a complex dimensionless number which can 
be expressed in terms of the gas thermal Peclet number 
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Peth = Ro^/Xg^ where Xg is the thermal diffusivity of the 



gas in ambient conditions [55l . |56 1 . 

Figure [T] displays the values calculated for 11* and n*j^ 
for an air bubble of ambient radius Rq = S /im driven 
at 20 kHz in water at ambient pressure and temperature. 
First, it is seen that the power dissipated either by viscous 
friction (thick solid line) or by thermal diffusion (thick 
dashed line) quickly rises in the neighborhood of the Blake 
threshold (where approximately 11* ~ IIJ'jj ~ 1), well above 
their value predicted by linear theory (between 5 and 6 
orders of magnitude). This clearly demonstrates the need 
for exact nonlinear bubble dynamics to calculate realistic 
values of the energy dissipated by inertial bubbles. 




Figure 1: Dimensionless power dissipated by an argon bubble of 
ambient radius Rq = 3 fim in water, at 20 kHz: by viscosity U* [thick 
solid line, from Eq. Ill3l l1: by thermal diffusion, n'j^ [thick dashed 
line, from Eq II12| |] . The thin lines are the corresponding values 
obtained from linear theory, Eqs. II15I I- II16II (solid; 11* jj^^; dashed: 
n*j^ lin)' vertical dash-dotted line represents the Blake threshold 
calculated by Eq. ^ . 



threshold with Blake threshold, Rozenberg's result sug- 
gests that nth+Hv would scale as P^ — P^ for P > Pb, and 
would be under the threshold. This is almost consistent 
with our results, except that redrawing Fig. [I] with linear 
scale (not shown) would reveal a linear dependence rather 
than a quadratic one. However, on one hand, Rozenberg's 
results apply to 500 kHz fields, and on the other hand, it is 
highly probable that the bubble density N also depends on 
the local sound field, despite we will consider N constant 
above the Blake threshold in the model developed below 
(see Sec. HI). 

We repeated the same calculation for a 8 bubble 
(Fig. [2]). The scale is chosen identical as Fig. [1] in order to 
make the comparison easier. The same conclusions apply 
except that the increase of viscous dissipation H* over the 
Blake threshold is lower than for the 3 /im bubble, and re- 
mains of the same order of magnitude as Hth for moderate 
driving pressures. 




Figure 2: Same as Fig. ^ for a 8 /im bubble. 



Another interesting feature is that, for the parameters 
used in Fig. [TJ viscous dissipation becomes much larger 
than the thermal one (more than one order of magnitude), 
for driving pressures above the Blake threshold, whereas 
linear theory predicts the opposite in this parameter range. 
Viscous dissipation in the liquid is thus found to largely 
predominate over the thermal one for 3 fim inertial bub- 
bles. 

It is also interesting to interpret these results in the 
light of the experimental data reported by Rozenberg [4l| , 
who fitted the volumic power dissipated in the cavitation 
zone by the following function of sound intensity /: 



V 



A{I- 
0, 



I>It 
I<It 



(17) 



where It is the intensity cavitation threshold. Identifying 
V with Af(nt,h + Hv), noting that sound intensity / scales 
as P^ for traveling waves, and identifying the cavitation 



To assess more clearly the dependence of Hv and Hth on 
the ambient radius Rq, we calculated Hv and Hth at con- 
stant P* = 1.5, but varying Rq. The result is displayed in 
Fig. 131 Viscous dissipation Hv is much larger than ther- 
mal dissipation Hth just above the Blake threshold, and 
decreases below Hth only above Rq ~ 10 ^m. 

More curves like the ones of Figs. [Tl [2] and |3] could be 
drawn, but we can summarize the comparison of Hv and 
Hth above the Blake threshold as follows: Hv predominates 
for larger drivings and smaller bubbles, while the opposite 
is true for larger bubbles and smaller drivings. Since for 
large drivings, surface instabilities maintain the ambient 
radii of inertial bubbles in a small interval just above the 
Blake threshold [H [s^l , this suggests that viscous fric- 
tion would be the predominant dissipation phenomenon in 
cavitation clouds. 

The real power dissipated by an inertial bubble is there- 
fore larger than the one predicted by linear theory by sev- 
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Figure 3: Same as Figs. [T] and [2] but varying for P* = 1.5. The 
vertical dasli-dotted line represents the Blake threshold. 



eral orders of magnitude. We therefore expect the real 
wave attenuation in a hquid containing inertial bubbles 
(above the Blake threshold) to be much higher than the 
value calculated by linear theory. We will quantify this 
point in Sec. [31 

Although the above results are sufhcient to carry on 
the development of our model, it is instructive to close 
this section by relating the dissipation functions Ilth and 
IIv to the conservation of acoustic energy, generalizing the 
conservation equation proposed in the original paper of 
Caflish et al. [32] ■ 

2.4. Conservation of energy in the bubbly liquid 

The term pdV/dt can be eliminated between equations 
(HJ and (O, by multiplicating the latter by N, to obtain a 
global energy conservation equation of the bubbly liquid: 



d_ 

di 



1 



1 2 



NKi 
dV 



AnNaR' 



(18) 



V. (pv) ^Npg—~ NlGnpiRR 



Equation (fTS|) represents the conservation of mechanical 
energy of the bubbly liquid; 

• p^/{2picf) is the elastic potential energy stored by 
the pure liquid involved in the propagation of the 
wave, 

• piv'^ /2 is the kinetic energy per unit volume of the 
pure liquid involved in the propagation of the wave, 

• NKi is the kinetic energy per unit volume of the 
liquid in its radial motion around the bubbles, 

• AnNaR^ is the interfacial potential energy per unit 
volume, 



• J3V is the acoustic intensity, or flux density of me- 
chanical energy. It is supplied at a vibrating bound- 
ary in contact with the bubbly liquid, typically by 
the oscillating motion of the sonotrode ^] . 

In what follows, we will assume periodic oscillations 
of all the fields. Averaging Eq. ((T5)) over one cycle, the 
time-derivative in the left side cancels and we get: 



V. (pv) = -N (Hth + Hv) . 



(19) 



Equation (|19p is the conservation of mechanical energy 
averaged over one period of oscillation, and has a clear 
physical interpretation: the balance between the acous- 
tic energy leaving a volume of bubbly liquid and the one 
reaching it is always negative, owing to thermal loss in the 
bubble and viscous friction in the radially moving liquid. 
Each bubble therefore appears as a dissipator of acoustic 
energy, owing to these two phenomena. The physical ori- 
gin of wave attenuation is thus self-contained in the Caflish 
model, even for nonlinear oscillations, provided that a cor- 
rect model is used to describe thermal diffusion in the 
bubble interior. Caflish and co-workers proposed a con- 
servation equation similar to (jl8p . disregarding viscosity 
and assuming isothermal oscillations, in which case me- 
chanical energy is conserved ^32^ . It should also be noted 
that Eq. (fTn]) reverts exactly the equation solved in ID by 
Rozenberg [41| in the case of purely traveling waves, but 
in the latter work, the dissipated power was fitted from 
experimental data, rather than being calculated ab initio 
from single bubble dynamics as done in the present work. 

3. The model 

3.1. Intuitive approach 

We first recall that the velocity field can be eliminated 
between Eq. ^ a nd (O ) to yield an equation involving only 
the pressure field [20|, [S^l : 



1 d^p d^P 
Pi- 



cf dt' di" 
Setting the pressure field p as a mono-harmonic wave: 
1 



(20) 



P(r,0 



P(r)e"^* + P{r)e- 



the linearization of the above equation and of the bubble 
dynamics equation allows to show that the complex field 
P fulfills an Helmholtz equation: 

V^P + k^P = 0, 

where the complex wave number is given by the linear 
dispersion relation [13, [l^, [H [H [H : 



n, — Q 



AttRqu^N 



2ibijj 



(21) 
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In Eq. (I2ip . wq is the resonance frequency and b the damp- 
ing factor, respectively given by 



2b = 



pRl 

Po(l + g) 



(22) 

(23) 



It can be readily seen, that, even for sub-resonant bubbles 
(w < Wo), the wave number is complex because of the 
damping factor 6, which, as expected from the discussion 
in Sec. 12.31 is correlated with the heat loss from the bubble 
and the viscous friction in the liquid. The imaginary part 
of the wave number represents the attenuation factor of the 
wave, and can be easily calculated by setting k — kr — ia 
and identifying kr and a from Eq. (j2ip . 

Generalizing this simple theory for inertial cavitation 
sounds unrealistic, since the bubble dynamics cannot be 
reasonably linearized for inertial oscillations. Thus all the 
fields are not mono-harmonic anymore and the problem 
cannot be reduced to an Helmholtz equation. However, 
for periodic oscillations, either linear or not, the correla- 
tion between the energy dissipated by each bubble over 
one cycle and the attenuation of the wave remains a uni- 
versal principle, formalized by Eq. and constitutes the 
guideline of the following derivation. 

We will therefore show that the first harmonic com- 
ponent of the field (at the frequency uj of the driving) 
approximately follows an Helmholtz equation, but whose 
wave number is directly expressed as functions of the dis- 
sipation functions Hth and Hv presented in the precedent 
section. This procedure allows to generalize the linear 
model, in the sense that the time-variable is eliminated, 
but keeping realistic values for the energy dissipated by 
inertial bubbles. 

3.2. Derivation of the model 

We decompose the pressure field into a sum of a time- 
average pressure Pm, a first harmonic pressure pi, oscillat- 
ing at the frequency of the ultrasonic source, and harmonic 
terms noted posc , that could be written as a Fourier series 
starting with a term at the frequency 2uj: 

p(r,i) ==_Pm(i") +pi{r,t) + 

Pose 

{r,t). (24) 
The first harmonic pressure field pi is expressed as: 



Pi(r,t) = i(P(r)e-*+P(r)e-^"* 



(25) 



where over-lines denote complex conjugate. Next, we set 
w the primitive of the first harmonic pressure field: 



w = l—[ P(r)e'"* - P(r)e" 
2 lu) ' 



(26) 



'^We assume here for simplicity that there is no subharmonics 
or ultra-harmonic terms, but the following reasoning can always be 
generalized by taking time-averages over the largest period of the 
pressure field. 



Multiplicating the propagation equation (1201) by w and 
averaging over one acoustic period yields: 



( WV P) = ^ { W 7T ) — Pl{ W 7T- 

^ ' cj\ dt^ / \ dt^ 



(27) 



Integrating by parts, using the definition of w, and the fact 
that all quantities are periodic, we obtain: 



(wW'^p) 



dp 
di 



'dt 



(28) 



Using the decomposition (1^^ . it can be easily seen that 
(w^'^p) = (luV^pi) and that {pidp/dt) = 0. Besides, 
using Eq. ((24| . the second term of the right-hand-side of 
((28|l can be expressed as: 



dp 



^'mrVTt 



dp 



dl3_ 

Pose ) , 



since {p^d/3/dt) = p^ {d/3/dt) — 0. We now make the em- 
pirical assumption that {poscdf3/dt) is negligible. A rigor- 
ous justification for this assumption is difficult in absence 
of results on the respective orders of magnitude of pi and 
Pose- However, unpublished measurements show that the 
latter is generally one order of magnitude lower than the 
former, so that for now, we assume that the assumption is 
justified. We therefore conclude that: 



Pi 



dt 



d/3 



(29) 



A physical interpretation of this approximate equation can 
be given by looking at Eqs. (H]), ([T]): it reverts to con- 
sider that the interaction between the acoustic field and 
the bubbles only occur through the first harmonic part of 
the field, and that the bubble mainly responds to this first 
harmonic content. We will term this hypothesis as "first 
harmonic approximation" (FHA). From this assumption 
and the above derivation, Eq. (|28p takes therefore the ap- 
proximate form: 



(wV^Pi) = piip 



dl_ 

'dt 



(30) 



and using Eq. ([SJ, Eq. ([50)1 . we finally obtain: 

(u;vV) = -A^(nth + nv). (31) 

We can now use the harmonic expressions and ([25)) of 
Pi and w, to obtain: 



4w 



(PV'P - PV^P) = -piN (Hth + 



and, dividing both sides of this equation by |Pp, P is 
finally found to fulfill: 



P 



(32) 
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where 9 denotes the imaginary part. We therefore see 
that if P were to fulfill an Helmholtz equation, the wave 
number would necessarily satisfy following relation: 

5(fc^)=-2p,.iV^^iii^. (33) 

This equation is a generalization of the linear case repre- 
sented by Eq. (PT|) . but here, Ilth and IIv can be estimated 
from fully nonlinear bubble dynamics. By the way, it can 
be checked after a few algebra that, linearizing Ilth and 
IIv, Eq. ([55)) yields the same results as taking the imagi- 
nary part of the dispersion relation (pij) . For linear oscilla- 
tions, Ilth and IIv scale as |Pp (see left part of the curves 
in Figs. [l]and[l]), so that linear theory yields a value of 
9 (fc^) independent of the driving amplitude This is 
no longer the case for nonlinear oscillations and Eq. 
yields a value of 5 (fc^) , which now depends on the local 
magnitude of the acoustic pressure \P\. 

The idea of the present paper is thus to use Eq. 
by using the nonlinear values of Ilth and IIv obtained in 
Sec. 12.31 to calculate 5(fc^), and, relying on Eq. to 
introduce the latter in a nonlinear Helmholtz equation: 

V2p + fc2(|p|)p = o. (34) 

Clearly, owing to the approximations made above, some 
additional terms would appear in the exact equation ful- 
filled by P. However, Eq. ([551) has the advantage to clearly 
link the attenuation factor to the real dissipation of en- 
ergy by the bubbles. Since it only yields the imaginary 
part of fc^, there remains the problem of calculating its 
real part. For now, we still use the linear dispersion re- 
lation to evaluate 9fi(fc^), and defer the discussion of this 
approximation below: 

cf OJ^ - UJ^ 

The attenuation coefficient and the real part of the wave 
number can now be deduced from: 

k = kr — ia, (36) 

and by identification with 

Figure |4] displays the attenuation coefficient a calcu- 
lated by following this procedure (thick solid line) , for 5 fim 
bubbles, and a typical 59] void fraction f^Q — 5 x I0~^, as 
a function of the acoustic pressure \P\. The attenuation 
coefficient rises abruptly for acoustic pressures just above 
the Blake threshold, as do Hth and Hv, and becomes about 
4 orders of magnitude larger than its linear value [thin 
solid line, calculated from Eq. (pij) ]. This demonstrates 
that a cloud of inertial cavitation bubbles damps out the 
incident wave much more drastically than linearly oscillat- 
ing bubbles. Moreover, contrarily to the linear prediction, 
the attenuation coefficient increases with the wave peak- 
amplitude. Thus, increasing the source vibration ampli- 
tude does not necessarily produce a more extended bubble 



field since increasing the acoustic pressure also increases 
the attenuation. This self-saturation phenomenon is well- 
known in cavitation experiments [36^ , and will be demon- 
strated in the simulations of the next section. 



1000 




'0 0.5 1 1.5 2 2.5 

p 

Figure 4: Real part (dashed) and imaginary part (solid) of the wave 
number k. The thin horizontal lines are predictions from linear the- 
ory II21I I and the thick lines are results calculated from Eqs. II33I I. I I35II . 
The vertical dash-dotted line represents the Blake threshold. 

The real part of the wave number is also displayed in 
Fig. m (thick dashed line), and the constant linear value 
predicted by (PT|) is recalled (thin dashed line) for com- 
parison. It is interesting to note that, above the Blake 
threshold, kr closely follows a. This comes from the fact 
that the imaginary part of A:^, calculated from the dissi- 
pation functions by Eq. is much larger in absolute 
value than its real part (j35p . This can be seen by writing 
the complex wave number as: 

k'^ = K^eiip[i{e-TT/2)], (37) 

where e is a small number, since 9 (A:^) is negative and 
large. The wave number k therefore reads: 

fc = Xexp[i(e/2-7r/4)], (38) 

and is therefore almost equal to K{1 — i)/^/2, so that we 
indeed have fcj. — a- 

The ratio a/kr has a strong physical sense. The atten- 
uation of the wave over one wavelength A is exp(— aA) = 
exp(— 27ra/fcr). Thus, if as in the present case a is of the 
same order of magnitude as kr , the attenuation of the wave 
over one wavelength is of the order of exp(— 27r) ~ 0.002. 
This means that as soon as the imaginary part of k^ is 
much larger than its real part, attenuation will play a dom- 
inant role whatever the precise value of its real value. This 
is why the precise choice of Sft (fc^) is of minor importance, 
and Eq. ([55)1 is a good compromise. 
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4. Results 

4.1. ID wave profiles 

We consider a tube of length L filled with water, bounded 
on the left by a piston which oscillating displacement reads: 



U{t) ~ Uq cosujt 



(39) 



and on the right by an infinitely soft boundary, imposing a 
zero acoustic pressure. This arbitrary boundary condition 
was chosen so that a standing wave should be obtained in 
the absence of bubbles. It can be easily changed to differ- 
ent and more complex conditions, as will be exemplified in 
the companion paper. 

We consider 5 /Ltm air bubbles. This choice is partially 
justified by experimental measurements of bubble size dis- 
tributions at low frequency 0, . In order to solve (IM]) 
along with Eqs. (|55)) - (|55|) . the bubble density N must be 
known. For now, we consider that bubbles are only present 
in the zones where the acoustic pressure amplitude is above 
the Blake threshold Eq. ([5]), and with a uniform density: 



N 



No ii\P\>PB 

if IPI < Pb 



(40) 



The nonlinear Helmholtz equation along with (I40p and 
the above boundary conditions is solved using the com- 
mercial COMSOL software, and a mesh convergence was 
performed. 

Figure [S] displays the profiles of the peak acoustic pres- 
sure \P*\ — P/po obtained for various amplitude of the 
source. For the smallest amplitude of the emitter Uq = 
0.2 /im, we recover a standing wave profiles in the pure 
liquid (dash-dotted line) . For a slightly larger vibration of 
the emitter Uq = 0.5 /im (dashed line), the acoustic pres- 
sure at the antinodes is just above the Blake threshold, so 
that the bubbles present here start to dissipate some en- 
ergy. This yields nonzero acoustic pressures at the nodes, 
but the profile remains globally similar to a linear standing 
wave profile. When the amplitude of the source is much 
larger (Uq = 5 /xm, solid line), the wave profile completely 
changes, and is drastically attenuated in a zone of about 
1 cm width near the emitter. This is due to the fact that 
the acoustic pressure near the emitter is larger than the 
Blake threshold, so that the bubbles present in this zone 
dissipate a lot energy. The remaining part of the profile is 
similar to a damped linear standing wave. 

In order to emphasize the importance of the nonlin- 
ear energy dissipation accounted for by our model, we 
present in Fig.lHla comparison of the upper profile of Fig. [S] 
(C/q = 5 fim, thick solid line), to the profile that would be 
obtained either by using the linear relation dispersion (PT|) 
with the same bubble density (thin solid line), or in the 
pure liquid (thin dashed line) . The important conclusion is 
that the two linear models predict unrealistic huge values 
of the acoustic pressure, while our model yields commonly 
measured amplitudes at 20 kHz (typically 1.5-3 bar [§|). 



|P*| 1 




0.02 0.04 0.06 0.08 

X (m) 



0.1 



Figure 5: Peak value of the dimensionless pressure field, calculated by 
solving numerically Eq. 1134 l l for various emitter displacement ampli- 
tudes. Solid line: Uo = 5 fim; dashed line: Uo = 0.5 /im; dash-dotted 
line: Uq = 0.2 fira. 



20 r 




0.04 0.06 
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Figure 6: Wave profiles for an amplitude of the emitter of 5 fim. 
Thick solid curve: predicted by the present model (same as the thick 
solid curve of Fig. [Sj; thin solid curve: obtained by the linear dis- 
persion relation Eq. II21I I: thin dashed curve: obtained in the pure 
liquid. 



4-. 2. Standing and traveling waves 

The phase 9 between the pressure field and the pressure 
gradient allows to determine whether the wave is travel- 
ing or standing. For a purely traveling wave (typically 



p{x,t) 



^i{ojt — kx) 



), pressure and pressure gradient are 



in phase quadrature, so that |sin0| = 1. Conversely, for 
a purely standing wave (typically p{x,t) ~ cos(/ca:)e'"*), 
pressure and pressure gradient are in phase or in phase op- 
position, so that I sin^l = in the latter case [60[. Thus, 
the quantity sin^ 9 can be used as a measurement of the 
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traveling character of the wave. 

In the configuration studied here, where the domain is 
closed with perfectly reflecting boundaries, linear acous- 
tics without dissipation would predict a perfect standing 
wave. However, if there is attenuation in the medium, a 
traveling wave component appears, because the reflected 
wave is of lower amplitude than the incident wave. This 
can be checked on Fig. [71 where sin^ 9 is displayed for the 
same simulation conditions as Fig. [SJ It is seen that for 
low driving amplitudes {Uq = 0.2 /im, dash-dotted line), 
sin^ 6* is everywhere, so that we have an almost perfect 
standing wave (which was clearly visible on Fig. [5]). But 
for higher emitter amplitude {Uq = 0.5 /im, dashed line), 
sin^ 9 starts to increase everywhere in the medium, espe- 
cially near the pressure antinodes, and for C/q = 5 /im (solid 
line) , sin^ 9 progressively increases toward 1 in a large part 
of the medium. 




0.02 0.04 0.06 0.08 0.1 

X (m) 

Figure 7: Phase between pressure and pressure gradient in the same 
conditions as Fig. [5] The hne-styles are the same as for Fig. [5] 

Finally, Fig.|S]confirms that, as shown above [see Eq. (1551) ]. 
the phase of the complex wave number k is close to — 7r/4 
in zones where the bubbles oscillate inertially. The wave 
number k is thus proportional to 1 — i, which means that 
the attenuation factor a and the real part A:,, of the wave 
number are of the same orders of magnitude. 

5. Conclusion 

Inertial bubbles dissipate much more energy than a lin- 
early oscillating bubble, both by thermal diffusion in the 
gas and viscous dissipation in the liquid, the latter mecha- 
nism being dominant for bubble ambient radii lower than 
10 /im. The wave attenuation in an inertial cavitation field 
is therefore much larger than the value predicted by the 
classical linear dispersion relation (by typically 4 orders 



O.OSr 




X (m) 

Figure 8: Phase of the complex wave number k divided by 7r, in the 
same conditions as Fig. (5] for Uq = 5 fim (sohd line) Uq = 0.5 ^m 
(dashed line), and Uo = 0.2 fim (dashed-dotted line). For the largest 
amplitude, the wave number near the emitter is seen to approach 
—it/4, as expected from Eq. i'dHt . 

of magnitude) . Although the latter conclusion is qualita- 
tively intuitive, to our knowledge, no quantitative estima- 
tion has ever been reported. 

Under the assumption that the bubbles are mainly ex- 
cited by the first harmonic content of the acoustic field, the 
latter fulfills approximately a nonlinear Helmholtz equa- 
tion. The imaginary part of the squared wave number is 
estimated rigorously from the energy dissipated by a sin- 
gle bubble, which can be easily calculated by solving a 
bubble dynamic equation. The real part is still arbitrarily 
estimated from the linear theory, but this arbitrary choice 
was shown to be of low importance, owing to the the huge 
value of the imaginary part. This has the importance con- 
sequence that in bubbly zones, the attenuation factor is of 
the same order of magnitude as the real part of the wave 
number, which results in a strong attenuation of the wave. 

The model has been solved in a typical ID-domain, 
and yields as expected a strongly attenuated wave pro- 
file near the emitter for high amplitude vibrations of the 
latter. The amplitude of the calculated acoustic pressure 
fields are realistic, contrarily to linear theory. This strong 
attenuation yields in turn a traveling component in the 
wave, where purely standing waves would be expected in 
a non dissipative medium enclosed by perfectly reflecting 
boundaries. 

It is interesting to note that following the present re- 
sults, attenuation, and therefore wave structures, are mainly 
governed by viscous dissipation involved in the bubble ra- 
dial motion, the thermal effects in the bubble playing a 
minor role. This conjecture might be checked experimen- 
tally by measuring the wave attenuation for solutions of 
different viscosities and with different dissolved gas. 

The choice of the incompressible Rayleigh-Plesset equa- 
tion to model the bubble dynamics may be questioned. 
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Although this is the original formulation of the Caflish 
model, the compressibility of the liquid produces sound 
scattering, and contributes therefore to the attenuation of 
waves in bubbly liquids, as is well-known in the linear case 
17, [2^. One may therefore replace Eq. ^ for example 
by a Keller equation 61-6^, and reformulate the energy 
equation (jS]) to exhibit an additional contribution of radi- 
ation Ha in its right-hand-side. This would in turn add a 
contribution in Eq. p3p , and produce more wave attenua- 
tion. However, the procedure is not straightforward, and 
the energetical interpretation in this case is less easy. One 
of the reasons for that is that compressible bubble dynam- 
ics equations are not exact solutions of the basic physical 
principles [g^, [l^ , but only first terms of expansions in the 
parameter R/q. It is also expected that sound scattering 
also modifies the real part of fc^, which again raises the 
issue of a correct expression for the latter. However, it 
may be conjectured that, in the low frequency range stud- 
ied here, the power loss by sound scattering is much lower 
than the one produced by viscous dissipation, because, as 
for thermal effects, sound radiation occurs mainly in the 
vicinity of the collapse. Thus we expect that the model in 
its present form catches the main dissipation phenomenon 
and that the values proposed for 3(fc^) is a good estima- 
tion. This will be examined in more details in future work. 

The occurrence of traveling waves, aside of the issue 
of the Bjerknes forces examined in the companion paper, 
may also have fundamental consequences on the final stage 
of the bubbles collapse. Indeed, it has been shown recently 
that bubbles in traveling waves are more exposed to shape 
instabilities and can undergo jetting, which reduces the fi- 
nal collapse temperature [6^, compared to a spherically 
collapsing bubble. This would therefore influence the esti- 
mation of the heat lost by a single bubble, but the spherical 
collapse model used in the present study yields an upper 
value. 

Besides, measurements of the acoustic field in conical 
structures has revealed the presence of a time-independent 
mean pressure field, which amplitude may be comparable 
with the first-harmonic part i36j. Our model does not 
catch this feature, and there is yet no correct theoreti- 
cal description of this phenomenon. We emphasize how- 
ever that our derivation of the imaginary part of the wave 
number is valid even in this case, since our decomposition 
of the field Eq. accounts a priori for the presence of 
such a mean field. This suggests that the present model 
could be supplemented by a specific equation describing 
this mean pressure field, which remains to be determined. 

To conclude, we believe that the present model opens 
the way to more realistic simulations of the coupled evo- 
lution of the cavitation field and the acoustic field. The 
nonlinear Helmholtz equation is relatively easy to solve 
and constitutes a viable solution halfway between a fully 
nonlinear simulation of the Caflish equations, which re- 
quires painful, if not intractable, temporal integration, and 
a fully linearized model which, as shown above, yields un- 
realistic acoustic pressure values. The companion paper 



will address the calculation of the Bjerknes forces in the 
acoustic fields calculated with the present model, and the 
resulting bubble structures predicted in more complex 2D 
configurations. 
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List of Figures 



1 Dimcnsionless power dissipated by an ar- 
gon bubble of ambient radius Rq — 5 fim 
in water, at 20 kHz: by viscosity IT* [thick 
solid line, from Eq. (IT51) ]: by thermal diffu- 
sion, n^j^ [thick dashed line, from Eq ([T^ ]. 
The thin lines are the corresponding values 
obtained from linear theory, Eqs. (fT5|) - (IT6|) 
(solid: n*j;jj; dashed: IIjjjjjjj). The ver- 
tical dash-dotted line represents the Blake 
threshold calculated by Eq. © 5 

2 Same as Fig. [T] for a 8 /xm bubble [5 

3 Same as Figs. [T] and [U but varying Rq for 
P* = 1.5. The vertical dash-dotted line rep- 
resents the Blake threshold 

4 Real part (dashed) and imaginary part (solid) 
of the wave number k. The thin horizontal 
lines are predictions from linear theory (I2ip 
and the thick lines are results calculated 
from Eqs. The vertical dash- 
dotted line represents the Blake threshold. . 

5 Peak value of the dimcnsionless pressure field, 
calculated by solving numerically Eq. 

for various emitter displacement amplitudes. 
Solid line: Uq — 5 ^m; dashed line: Uq = 
0.5 ^m; dash-dotted line: Uq = 0.2 /xm. . . 

6 Wave profiles for an amplitude of the emit- 
ter of 5 /xm. Thick solid curve: predicted by 
the present model (same as the thick solid 
curve of Fig. [5]); thin solid curve: obtained 
by the linear dispersion relation Eq. (|2ip : 
thin dashed curve: obtained in the pure liq- 
uid i 

7 Phase between pressure and pressure gradi- 
ent in the same conditions as Fig. [5j The 
line-styles are the same as for Fig. [5] . . . . [13 

8 Phase of the complex wave number k di- 
vided by 77, in the same conditions as Fig.[51 
for Uq = 5 fiiR (solid line) Uq = 0.5 /zm 
(dashed line), and Uq = 0.2 (dashed- 
dotted line). For the largest amplitude, the 
wave number near the emitter is seen to ap- 
proach — 7r/4, as expected from Eq. 
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